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It is demonstrated how to generate time series with tailored nonlinearities by inducing well-defined 
constraints on the Fourier phases. Correlations between the phase information of adjacent phases 
and (static and dynamic) measures of nonlinearities are established and their origin is explained. By 
applying a set of simple constraints on the phases of an originally linear and uncorrelated Gaussian 
time series, the observed scaling behavior of the intensity distribution of empirical time series can be 
reproduced. The power law character of the intensity distributions being typical for e.g. turbulence 
and financial data can thus be explained in terms of phase correlations. 
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Introduction. The clearest yet most general definition 
of nonlinearity in time series g{t) is given in the Fourier 
representation 

1 

G{k) = FT{g{t)) = - ^ (1) 

t=0 

of the data. Linear time series are fully characterized by 
the modulus \G{k)\ of the complex valued Fourier coef¬ 
ficients G{k) = \G{k)\e'^'^^^\ while the phases (l){k) are 
uncorrelated and uniformly distributed in the interval 
(p G [—7r;7r]. Any nonlinearity is coded in the Fourier 
phases (l){k) and correlations among them. Deviation 
from the randomness of the phases are thus equivalent 
to the presence of nonlinearities in the time series. As 
yet, only little attention has been paid so far to the ex¬ 
plicit analysis of the information contained in the Fourier 
phases to characterize nonlinearities, although a lot of in¬ 
sights about nonlinearities may be gained by better un¬ 
derstanding the meaning of phases. 

The definition of nonlinearity via the randomness of 
Fourier phases is - on the other hand - at the heart of 
algorithms for generating so-called surrogate data sets, 
which were developed to test for weak nonlinearities in 
a model-independent way [T]. These surrogates are sup¬ 
posed to have the same linear properties as a given data 
set, while all nonlinearities are wiped out. The removal 
of the nonlinear correlations is achieved by replacing the 
phases 0(/c) with a set of uncorrelated and uniformly dis¬ 
tributed ones. Refinements of the Fourier-based methods 
for generating surrogates aimed at preserving both the 
power spectrum \G{k)\‘^ and the amplitude distribution 
of the time series g{t) in real space [IHil- The addition of 
(iterative) rank-ordered remapping of the phase random¬ 
ized data onto the original amplitude distribution led to 
surrogates with the desired amplitude spectrum mil]. 
Applying the lAAFT method in the wavelet domain al¬ 
lowed for the generation of surrogates which also preserve 
the local mean and variance of the original signal mm. 
However, it was found recently that these (iterated) am¬ 
plitude adjusted ( (I) A AFT ) surrogates may not be lin¬ 
ear, since the randomness of the phases is guaranteed 


only before the first remapping step. One can rather find 
phase correlations in surrogate realizations that may re¬ 
sult in a non-detection of nonlinearities in time series [5] . 
But this obvious flaw of (I) A AFT surrogates became a 
virtue as significant correlations between phase statistics 
and a measure for nonlinearity were found for the first 
time (see [5] and insets in Fig. . 

Connections between correlations among Fourier phases 
and higher order statistics could also identified by analyz¬ 
ing the cosmic microwave background radiation (CMB). 
Several studies of both the WMAP and PLANCK data 
involving surrogates revealed that there are phase cor¬ 
relations at large scales in the CMB which lead to pro¬ 
nounced anisotropies (see e.g. IHHI). Recently it was 
demonstrated that the observed phase correlations can 
gradually be diminished when subtracting suitable best- 
fit (Bianchi-)template maps. The weaker phase correla¬ 
tions lead in turn to a vanishing signature of anisotropy 
as identified with higher order statistics [9]. 

The relations between phase information and higher or¬ 
der statistics in (I) A AFT surrogates and the CMB data 
were only found in a heuristic manner. 
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FIG. I: Gaussian random uncorrelated noise (black). The 
colored points show time series with linear phase correlations 
among adjacent phases (A = I) with cr^ = I.O and = 0 
(blue) and = tt (red). 

To allow for a systematic investigation of phase cor- 
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FIG. 2: Upper Row: Nonlinear prediction error '0 (upper left) and average connectivity (upper right) versus phase correlation 
coefficient c(A = 1) for 900 time series with imposed phase correlations derived from Gaussian noise as input time series g{t). 
Lower row: Same as upper row but with the Mrk 766 X-ray observation as input time series g{t). The data from revolution 
999 was binned with a binsize of 50 s leading to a time series with 1540 points. The insets show the corresponding results for 
400 lAAFT surrogates. 


relations and their corresponding nonlinearities in time 
series, it is desirable to start with the phases, constrain 
their correlations in a tunable and reproducible way and 
study the effects on the nonlinear statistics. 

Here, we present a method to generate time series with 
such tailored nonlinearities by imposing well-defined cor¬ 
relations on the Fourier phases and demonstrate how de¬ 
viations from linearity can be understood in terms of 
phase information. 

Methods. To address the relationship between phase cor¬ 
relations and measures for nonlinearity we calculate the 
nonlinear prediction error (NPLE) [10] as an example 
for a dynamical complexity measure with a good overall 
performance m and the average connectivity of (recur¬ 
rence) networks as an example for a structural complex¬ 
ity measure nauB]. The calculation of both measures 
relies on the representation of the time series in an ar¬ 
tificial phase space, which is obtained using the method 
of delay coordinates M- This is accomplished by us¬ 
ing time delayed versions of the observed time series as 
coordinates for the embedding space. The multivari¬ 
ate vectors in the d-dimensional space are expressed by 
9t = {9t,9t+T,gt+2T, ■ ■ .,9t+{d-i)T) , where r is the delay 
time and gt denotes the value of the (discretized) time 
series at time step t. 

The comparison of the predicted behavior of the embed¬ 
ded time series based on the local neighbors with the 


real trajectory of the system leads to the definition of 
the NLPE -0 as 


^ = ij{d,T,T,N) 


( 2 ) 




M-l-T 


^ [g,+T-F{g,)]\ 


where F is a locally constant predictor, M is the length 
of the time series, and T is the lead time. The predic¬ 
tor F is calculated by averaging over future values of the 
= d + 1 nearest neighbors in the delay coordinate rep¬ 
resentation. We found that -0 remains rather constant for 
T>5, thus a value of T = 5 was used for this study. The 
dimension of the embedding space d and the delay time r 
have to be set appropriately. Since the time series of an 
Active Galactic Nuclei (AGN) being studied in the fol¬ 
lowing consists of less than 1600 data points, we use a low 
embedding dimension d = 3. Due to the long correlation 
time of this time series, we chose a relatively large delay 
time r = 250 according to the criterion of zero crossing 
of the autocorrelation function m- To allow for a direct 
comparison, we use the same values d = 3 and r = 250 
for the other time series with imposed phase correlations. 
The structural complexity of a time series with a limited 
number of points can be characterized with recurrence 
networks m- They are based on recurrence plots dl, 
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which describe how often pairs of points of a time series 
in the embedding space representation come close to each 
other. Linking such nearby points in a network represen¬ 
tation of the data and omitting self-loops leads to the 
definition of the adjacency matrix Ai^j of the recurrence 
network [T^ 


Aij = e(e-||gi-gjil)(3) 

where the {gn} are the data points in embedding space 
and e is an appropriate threshold. Ai^j contains the whole 
information about the network. A common measure for 
the topological structure of the network is the average 
connectivity n which is calculated by 

^ n—1 

(“) 

^ ^ zy=0 


where ky = is the degree of node u. If the 

attractor of the nonlinear system is reconstructed with 
appropriate embedding parameters this network measure 
can be used as a test for nonlinearity. The threshold e 
is chosen such that k = 0.01 for the original time series. 
The same threshold is then used for the Gaussian time 
series with imposed phase correlations. 

To get a visual impression of correlations among the 
Fourier phases it is convenient to make use of so-called 
phase maps [18]. A phase map is defined as a two- 
dimensional set of points G = { 0 /c 5 0/c+a} where (j)^ is 
the phase of the mode of the Fourier transform and 
A a frequency delay. To quantify the degree of correla¬ 
tion between the phases 0 and 0 + A we calculate the 
correlation coefficient c(A), 


c(A) = 


{(p{k)4>{k + A)) 

^0(/e)^0(/c+A) 


( 5 ) 


Note that by using c(A) as correlation measure we re¬ 
strict ourselves to the simplest way of quantifying corre¬ 
lations among the phases that is only sensitive to linear 
correlations. 

Time series with phase eorrelations. As outlined in [7j , 
(I) A AFT surrogates can contain phase correlations lead¬ 
ing to statistically significant high or low values of c(A). 
A closer look at the corresponding phase maps reveals 
that the (anti-)correlations originated from stripe-like, 
patterns along the diagonal (i.e. with slope of one) or 
shifted relative to it. These patters thus indicate that 
phase pairs are linearly correlated with each other. One 
can further notice that for the time series stemming from 
X-ray observations of the AGN Mrk 766 the phase cor¬ 
relations are most pronounced for A = 1 . We reproduce 
such signatures by imposing correlations in the phase dis¬ 
tribution in the following way: The values for the phases 
0 (/c) are iteratively determined by relating A) with 
Hk) by 


(6) 


with d4> being a shift constant ranging from — tt to tt and 
T] describing a (Gaussian) noise term with given stan¬ 
dard deviation cr^. In the phase map picture 77 controls 
the width of the stripes and dcj) defines its position. The 
iteration is performed over the frequencies /c, where kg 
denotes the starting value and dk the step size of the 
iteration. 0 (/Cs) is drawn from a uniform distribution 
within the interval [— 7 r, 7 r]. The same is true for 0(/c) if 
this phase has not been set in a previous iteration step. 
In our first example we are interested in only correlating 
adjacent phases. Thus we apply Eq. j^with A = 1 to a 
Gaussian time series with zero mean and standard devi¬ 
ation of one. The step size is chosen to dk = 2. Thereby 
every phase is correlated to exactly one other phase for 
A = 1 , while the phases are not correlated for any fre¬ 
quency delay A greater than one. 

Fig.[T] shows how these phase correlations alter the time 
series. It becomes clearly visible that the correlations 
of adjacent phases induce fluctuations of the variance. 
Specifically, one recognizes a time interval where the fluc¬ 
tuations are larger than for the noise and another region 
where the fluctuations are smaller. Note that the over¬ 
all mean and standard deviation of the time series are 
exactly preserved since the power spectrum is kept con¬ 
stant. The shift constant controls the position of the 
region with higher fluctuations. If = 0, this region is 
located in the middle of the time series and it shifts to¬ 
wards the ends of the time series when dcj) approaches Ttt. 
By testing different values of A we further found that the 
number of regions with high fluctuations is given by the 
value for A. In Fig. we show the nonlinear predic¬ 
tion error -0 and the mean connectivity as a function 
of the correlation coefficient c(A = 1 ). The results are 
displayed for time series with imposed phase correlations 
(only) for A = 1 and varying d0 as derived from Gaus¬ 
sian noise and from the X-ray observation of the AGN 
Mrk 766. One can see that the shift constant d4> con¬ 
trols the (anti-)correlations of the phases. More impor¬ 
tantly, it becomes obvious that both the nonlinear pre¬ 
diction error and the mean connectivity are highly (anti- 
) correlated with the phase correlations as measured with 
c(A). Knowing that dcj) also controls the position of the 
regions with higher and lower fluctuations, we can now 
get a much more detailed understanding of how the phase 
correlations influence the calculation of the NLPE and 
the average connectivity. The embedding with the delay 
time of r = 250 in three dimensions leads to a truncation 
of the last part of the time series. Depending on whether 
the remaining time series has lower {dcj) ~ Ttt) or larger 
{dcj) ~ 0 ) fluctuations, one obtains larger or lower values 
for the NLPE leading to the observed anti-correlation be¬ 
tween psi and c. Similarly, lower fluctuations in the time 
series lead to a more connected recurrence network and 
vice versa, correlating n and c(A). 

In a second example we extend the formalism to gener¬ 
ate nonlinear time series with well-defined nonlinearities 
by simultaneously imposing linear phase correlations for 
a set of different frequency delays A. This is achieved 


0(/c + A) = 0(/c) -h # + 7 / 
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TABLE I: Parameters defining the imposed phase correlations 



A 

dk 

ks 

dcj) 

9 

1. Iteration 

1 

2 

1 

3.1415 

0.1 

2. Iteration 

3 

3 

1 

3.0 

0.08 

3. Iteration 

3 

3 

2 

3.0 

0.3 

4. Iteration 

5 

5 

2 

1.4 

0.2 

5. Iteration 

7 

7 

3 

3.1415 

0.25 

6. Iteration 

50 

50 

2 

3.0 

0.1 
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FIG. 3: Time series g{t) that is obtained from white Gaussian 
noise by imposing a set of six different linear phase correla¬ 
tions. 


by iteratively applying Eq. starting with low values 
of A and then proceeding to higher ones. Iterating over 
increasing frequency delays A ensures that phase cor¬ 
relations that were imposed in previous iteration steps 
are at least in parts preserved when new constraints for 
phase correlations at larger A are added. Table [l| sum¬ 
marizes the parameters for the six iterations used in our 
example. Fig. shows the time series which is obtained 
when the six constraints on phase correlations are im¬ 
posed on white Gaussian noise. One has to note that 
the time series has no linear correlation as the modulus 
\G{k)\ of the Fourier transform of the original random 
time series is left untouched. In the time series with 
phase correlations one can clearly identify a number of 
time intervals with larger fluctuations whose number, po¬ 
sition and strength are controlled by the parameters A, 
dk and 77 , respectively. The statistical properties of such 
a nonlinear time series can thus be tailored in a refined 
manner. The time series in this example was generated 
such that it resembles data often observed in economic 
time series m, where especially data from stock indices 
show intermittent behavior, i.e. extreme events, patterns 
of volatility clustering and phase correlations, while the 
autocorrelation vanishes. The distribution of the fluctu¬ 
ations is further analyzed by calculating the cumulative 
probability distribution P{g) of the normalized positive 
and negative values of g{t) (see Fig. [^. We find the 
expected leptokurtic distribution whose tail can be fit¬ 
ted with a power law P{g) ^ g~^ with a = 3.25 ± 0.16 



FIG. 4: Gumulative probability distribution P{g) of the nor¬ 
malized positive (black) and negative (red) values of g(t). The 
black dashed and red dotted lines show the respective distri¬ 
butions for the initial white Gaussian noise. 



FIG. 5: Gumulative probability distribution P(v) of the nor¬ 
malized volatility v. The black dashed line shows the respec¬ 
tive distributions for the initial white Gaussian noise. 


for the positive tail and a = 3.37 ± 0.47 for the negative 
tail in the region 5 < ^ < 20. These numbers are in re¬ 
markable agreement with those obtained from empirical 
studies of market indices [20]. We further studied the 
statistical properties of the volatility v{t) as defined as 
the average of \g{t)\ over a time window of length A, i.e. 
v{t) = 1/A Fig- shows the cumulative 

probability distribution P{v) for A = 5. As expected 
we find a distribution with fat tails, which can be fit¬ 
ted by P{v) v~^ with [3 = 2.93 ± 0.49 in the region 
5 < 7 ; < 20. Again, this is in very good agreement with 
the scaling properties of the volatility of price fluctua¬ 
tions observed in empirical data m- 
Finally, we note that analogies between price dynam¬ 
ics of market indices and the velocity differences in 
three-dimensional fully developed turbulence have been 
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pointed out by several authors (see e.g. [22l[23]). Conse¬ 
quently, the fat tails in the probability density functions 
of turbulence data may be also be understood in terms of 
phase correlations allowing for a better characterization 
and discrimination of different scenarios of turbulence. 

Summary. We have presented a new method to gen¬ 
erate time series with well-defined nonlinearities by im¬ 
posing linear correlations among the Fourier phases. We 
have shown that the phase correlations between adja¬ 
cent phases are tightly related with higher order statis¬ 
tics being estimated for the time series. These ”Wiener- 
Khinchin-like” connections between phase information 
and higher order statistics are to a large extent indepen¬ 


dent of the input time series. Furthermore, the scaling 
of fluctuation and of the volatility of a time series can 
be understood in terms of a set of linear phase correla¬ 
tions. We expect that further studies with time series 
and also spatial structures with tailored nonlinearities, 
for which not only linear but more complex constraints 
on the Fourier phases are imposed, will shed more light 
on both the meaning of Fourier phases and the different 
kinds of nonlinearities as they are observed in nature. 
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